Trophoblasts_hypoxia_pseudotime_fig2
Jose J Moreno-Villena
2025-04-09
Knitted doing: rmarkdown::render(“~/Documents/ucl_postdoc/FMI_supplementary_information/Scripts_to_upload_to_github/Planceta_STB_pseudotime_hypoxia_9April2025.Rmd”)
library(dplyr)
library(slingshot)
library(patchwork)
library(slingshot)
library(ggrepel)
library(condiments)
library(tidyr)
library(RColorBrewer)
library(Seurat)Subset Seurat object
all_fmi<-readRDS(file = "/Users/jose/Documents/ucl_postdoc/trophoblasts/FMI-all-20patients-20240827.rds")
tropho <- subset(all_fmi, CellTypeManual.l1 == "Trophoblast")
rm(all_fmi)
nojj<-unique(tropho$sample_id)[!grepl("F2044JJ",unique(tropho$sample_id))]
tropho <- subset(tropho, sample_id %in% nojj)
# subset to placenta
placenta_fmi<-subset(tropho, Tissue == "PVBP")
rm(tropho)
tropho_stb <- subset(placenta_fmi, CellTypeManual.l3 %in% c("STB","STB-precursor","CTB","CTB-proliferating"))
rm(placenta_fmi)Generate UMAP
tropho_stb<- FindVariableFeatures(tropho_stb, selection.method = "vst")
tropho_stb <- ScaleData(tropho_stb)## Centering and scaling data matrixtropho_stb <- RunPCA(tropho_stb, verbose = FALSE)
# UMAP calculation
tropho_stb <- RunUMAP(tropho_stb, dims = 1:10)## 17:12:21 UMAP embedding parameters a = 0.9922 b = 1.112## 17:12:21 Read 21465 rows and found 10 numeric columns## 17:12:21 Using Annoy for neighbor search, n_neighbors = 30## 17:12:21 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 17:12:23 Writing NN index file to temp file /var/folders/t2/8wy_rh2s7lv64wdz4bbh47km0000gp/T//RtmpdZwNq8/file17a9d1502fa53
## 17:12:23 Searching Annoy index using 1 thread, search_k = 3000
## 17:12:27 Annoy recall = 100%
## 17:12:28 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 17:12:29 Initializing from normalized Laplacian + noise (using RSpectra)
## 17:12:29 Commencing optimization for 200 epochs, with 881908 positive edges
## 17:12:35 Optimization finished# Create the DimPlot with custom colors
p_celltype<-DimPlot(tropho_stb, reduction = "umap", group.by = "CellTypeManual.l3", label = TRUE, label.size = 3, repel = TRUE, split.by = "Condition")
p_celltype
#####
#convert to SCE
seu_placenta_initial <- as.SingleCellExperiment(tropho_stb, assay = "RNA")Infer trajectory & pseudotime
seu_placenta <- slingshot(seu_placenta_initial, reducedDim = 'UMAP',
clusterLabels = colData(seu_placenta_initial)$CellTypeManual.l3,start.clus = "CTB-proliferating")Convert to data.frame to be used in ggplot
seu_placenta1 <- bind_cols(
as.data.frame(reducedDims(seu_placenta)$UMAP),
as.data.frame(colData(seu_placenta)[1:40]),
## need to remove all the slingshot and and tradeseq info. tables and can't be conversted into datafram "fitgam" and "crv" fields
slingPseudotime(seu_placenta)
) %>%
sample_frac(1)
# Create the initial ggplot
p <- ggplot(seu_placenta1, aes(x = umap_1, y = umap_2)) +
geom_point(aes(fill = Condition), col = "grey70", shape = 21) +
theme_minimal() +
scale_fill_viridis_d()
# Add the geom_path layer to the ggplot
curves <- slingCurves(seu_placenta, as.df = TRUE)
p1 <- p +
geom_path(data = curves %>% arrange(Order),
aes(group = Lineage, col = as.character(Lineage)), size = 1.5, arrow = arrow(type = "closed", length = unit(0.15, "inches"))) +
#theme_minimal() +
labs(col = "Lineage")
p1
Imbalance scores in condiments
scores <- condiments::imbalance_score(
Object = seu_placenta1 %>% select(umap_1, umap_2) %>% as.matrix(),
conditions = seu_placenta1$Condition,
k = 20, smooth = 40)
seu_placenta1$scores <- scores$scores
seu_placenta1$scaled_scores <- scores$scaled_scoresThe local scores are quite noisy so we can then use local smoothers
to smooth the scores of individual cells. The smoothness is dictated by
the smooth argument. Those smoothed scores were also
computed using the imbalance_score function.
# Create the initial ggplot
p <- ggplot(seu_placenta1, aes(x = umap_1, y = umap_2, col = scaled_scores)) +
geom_point() +
scale_color_viridis_c(option = "C")
# Add the geom_path layer to the ggplot
#curves <- slingCurves(seu_placenta, as.df = TRUE)
p1 <- p +
geom_path(data = curves %>% arrange(Order),
aes(group = Lineage), col = "green", size = 1.5,
arrow = arrow(type = "closed", length = unit(0.15, "inches"))) +
theme_classic() +
labs(col = "Imbalance") +
geom_point(data = curves %>%
group_by(Lineage) %>%
slice(round(seq(1, n(), length.out = 9))),
aes( size = 4),col = "green") +
theme_classic() +
labs(col = "Imbalance")
p1
#ppp4<-ggplot(seu_placenta1, aes(x = umap_1, y = umap_2, col = scaled_scores)) +
# geom_point() +
# scale_color_viridis_c(option = "C")
#ppp4Imbalance scores in Early
seu_placenta1_early<-subset(seu_placenta1, GA_Category %in% "Early")
scores <- condiments::imbalance_score(
Object = seu_placenta1_early%>% select(umap_1, umap_2) %>% as.matrix(),
conditions = seu_placenta1_early$Condition,
k = 20, smooth = 40)
seu_placenta1_early$scores <- scores$scores
seu_placenta1_early$scaled_scores <- scores$scaled_scores
# Create the initial ggplot
p <- ggplot(seu_placenta1_early, aes(x = umap_1, y = umap_2, col = scaled_scores)) +
geom_point() +
scale_color_viridis_c(option = "C")
# Add the geom_path layer to the ggplot
curves <- slingCurves(seu_placenta, as.df = TRUE)
p1 <- p +
geom_path(data = curves %>% arrange(Order),
aes(group = Lineage), col = "green", size = 1.5,
arrow = arrow(type = "closed", length = unit(0.15, "inches"))) +
theme_classic() +
labs(col = "Imbalance") +
#scale_color_brewer(palette = "Set1")+
geom_point(data = curves %>%
group_by(Lineage) %>%
slice(round(seq(1, n(), length.out = 9))),
aes( size = 4),col = "green") +
theme_classic() +
labs(col = "Imbalance")
#scale_color_brewer(palette = "Set1")
p1
#ppp4<-ggplot(seu_placenta1, aes(x = umap_1, y = umap_2, col = scaled_scores)) +
# geom_point() +
# scale_color_viridis_c(option = "C")
#ppp4Imbalance scores in Late
seu_placenta1_late<-subset(seu_placenta1, GA_Category %in% "Late")
scores <- condiments::imbalance_score(
Object = seu_placenta1_late%>% select(umap_1, umap_2) %>% as.matrix(),
conditions = seu_placenta1_late$Condition,
k = 20, smooth = 40)
seu_placenta1_late$scores <- scores$scores
seu_placenta1_late$scaled_scores <- scores$scaled_scores
# Create the initial ggplot
p <- ggplot(seu_placenta1_late, aes(x = umap_1, y = umap_2, col = scaled_scores)) +
geom_point() +
scale_color_viridis_c(option = "C")
# Add the geom_path layer to the ggplot
curves <- slingCurves(seu_placenta, as.df = TRUE)
p1 <- p +
geom_path(data = curves %>% arrange(Order),
aes(group = Lineage), col = "green", size = 1.5,
arrow = arrow(type = "closed", length = unit(0.15, "inches"))) +
theme_classic() +
labs(col = "Lineage") +
#scale_color_brewer(palette = "Set1")+
geom_point(data = curves %>%
group_by(Lineage) %>%
slice(round(seq(1, n(), length.out = 9))),
aes( size = 4),col = "green") +
theme_classic() +
labs(col = "Lineage")
#scale_color_brewer(palette = "Set1")
p1
#ppp4<-ggplot(seu_placenta1, aes(x = umap_1, y = umap_2, col = scaled_scores)) +
# geom_point() +
# scale_color_viridis_c(option = "C")
#ppp4Combine imbalance scores Early and Late
seu_placenta1_GA<-rbind(seu_placenta1_late, seu_placenta1_early)
# Create the initial ggplot
p <- ggplot(seu_placenta1_GA, aes(x = umap_1, y = umap_2, col = scaled_scores)) +
geom_point() +
scale_color_viridis_c(option = "C") + facet_wrap(~GA_Category)
# Add the geom_path layer to the ggplot
curves <- slingCurves(seu_placenta, as.df = TRUE)
p1 <- p +
geom_path(data = curves %>% arrange(Order),
aes(group = Lineage), col = "green", size = 1.5,
arrow = arrow(type = "closed", length = unit(0.15, "inches"))) +
theme_classic() +
labs(col = "Lineage") +
#scale_color_brewer(palette = "Set1")+
geom_point(data = curves %>%
group_by(Lineage) %>%
slice(round(seq(1, n(), length.out = 9))),
aes( size = 4),col = "green") +
theme_classic() +
labs(col = "Lineage")
#scale_color_brewer(palette = "Set1")
p1
#ppp4<-ggplot(seu_placenta1, aes(x = umap_1, y = umap_2, col = scaled_scores)) +
# geom_point() +
# scale_color_viridis_c(option = "C")
#ppp4seu_placenta1_GA<-rbind(seu_placenta1_late, seu_placenta1_early)
seu_placenta1_GA_nocbt<-subset(seu_placenta1_GA, umap_1 < 5)
# Create the initial ggplot
p <- ggplot(seu_placenta1_GA_nocbt, aes(x = umap_1, y = umap_2, col = scaled_scores)) +
geom_point() +
scale_color_viridis_c(option = "C") #+ facet_wrap(~GA_Category)
# Add the geom_path layer to the ggplot
curves <- slingCurves(seu_placenta, as.df = TRUE)
curves2<-subset(curves, Order > 38)
p1 <- p +
geom_path(data = curves2 %>% arrange(Order),
aes(group = Lineage), col = "green", size = 1.5,
arrow = arrow(type = "closed", length = unit(0.15, "inches"))) +
theme_classic() +
labs(col = "Lineage") +
#scale_color_brewer(palette = "Set1")+
geom_point(data = curves2 %>%
group_by(Lineage) %>%
slice(round(seq(1, n(), length.out = 7))),
aes( size = 4),col = "green") +
theme_classic() +
labs(col = "Lineage")
#scale_color_brewer(palette = "Set1")
p1
Plot pseudotime per cell in UMAP
# Create the initial ggplot
p2 <- ggplot(seu_placenta1, aes(x = umap_1, y = umap_2)) +
geom_point(aes(fill = Lineage1), col = "grey70", shape = 21) +
#theme_minimal() +
scale_fill_viridis_c()
# Add the geom_path layer to the ggplot
p3 <- p2 +
geom_path(data = curves %>% arrange(Order),
aes(group = Lineage, col = as.character(Lineage)), size = 1.5, arrow = arrow(type = "closed", length = unit(0.15, "inches"))) +
#theme_minimal() +
labs(col = "Lineage") +
scale_color_brewer(palette = "Set1")+
geom_point(data = curves %>%
group_by(Lineage) %>%
slice(round(seq(1, n(), length.out = 9))),
aes(col = as.character(Lineage)), size = 4) +
theme_classic() +
labs(col = "Lineage") +
scale_color_brewer(palette = "Set1")## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.# Print the updated plot
print(p3)
Function to calculate differential abundance in pseudotime
DAbundance_sum_expression_in_pseudotime <- function(df_sample, lineages = "Lineage1", gest_age = NULL, conditions = c("Control", "PE"), pseudotime_intervals = 1, tropho,genestoplot,expr_factor = 100) {
# Subset data frame by lineages and gestational age if specified
if (!is.null(gest_age)) {
### make the script:
#lineages = "Lineage3"
#gest_age = "Early"
#pseudotime_intervals = 1
#conditions = c("Control", "PE")
dfl1 <- df_sample[df_sample$lineages == lineages & df_sample$GA_Category == gest_age, ]
} else {
dfl1 <- df_sample[df_sample$lineages == lineages, ]
}
# Define the range of pseudotime values
min_pt <- min(dfl1$pseudotime)
max_pt <- max(dfl1$pseudotime)
# Create a sequence of pseudotime values at intervals of specified units
pseudotime_seq <- seq(min_pt, max_pt, by = pseudotime_intervals)
# Initialize density_df outside the loop if not already done
# Initialize empty dataframe to store density values for control and PE samples
density_df <- data.frame(Pseudotime = pseudotime_seq,
n_samples_Control = rep(0, length(pseudotime_seq)),
n_samples_PE = rep(0, length(pseudotime_seq)),
Density_Control = rep(0, length(pseudotime_seq)),
Density_PE = rep(0, length(pseudotime_seq)),
Wilcox_Pvalue = rep(0, length(pseudotime_seq)),
Mean_Control = rep(0, length(pseudotime_seq)),
Mean_PE = rep(0, length(pseudotime_seq)),
expr_score = NA,
done = "no"
)
# Loop through pseudotime sequence
for (i in seq_along(pseudotime_seq)) {
# Subset data for each pseudotime interval
subset_data <- dfl1[dfl1$pseudotime >= pseudotime_seq[i] & dfl1$pseudotime < pseudotime_seq[i] + pseudotime_intervals, ]
# Check if there are any cells in this interval
if (nrow(subset_data) == 0 || length(intersect(rownames(subset_data), colnames(tropho))) == 0) {
message(paste("No cells found for interval", i, ". Skipping..."))
density_df$done[i] <- "yes"
density_df[i,]<-NA
next # Skip to the next iteration
}
# Subset the Seurat object
interval_seurat <- subset(tropho, cells = intersect(rownames(subset_data), colnames(tropho)))
# Check if the subset Seurat object has any cells
if (ncol(interval_seurat) <= 1) {
message(paste("No cells remained after subsetting for interval", i, ". Skipping..."))
next # Skip to the next iteration
}
#genes_of_interest <- list(c("LEP", "LTF", "DERL3"))
#interval_seurat <- AddModuleScore(interval_seurat,
# features = genes_of_interest,
# name = "GeneScore")
#density_df$expr_score[i] <- mean(interval_seurat$GeneScore1)
# Get the expression data for genes of interest
#genes_of_interest <- c("HTRA1","FSTL3","LTF","LEP","FLT1")
genes_of_interest <- genestoplot
gene_expr <- GetAssayData(interval_seurat, slot = "data")[genes_of_interest, ]
# Calculate the sum of expression of all the genes for each cell
if(length(genes_of_interest) == 1 ) {
sum_expr <- t(as.data.frame(gene_expr))} else {
sum_expr <- colSums(gene_expr)}
# Calculate the mean of the sum across all cells in this interval
density_df$expr_score[i] <- mean(sum_expr)
density_df$done[i] <- "yes"
#}
# ggplot(data = density_df, aes(x = pseudotime)) +
# geom_line(aes(y = 10*expr_score, color = 'coral'), size = 1) +
# labs(x = "Pseudotime", y = "Mean Hypoxia score", title = paste(lineages,"and",gest_age)) +
# theme_minimal() +
# theme(legend.position = "top") +
# guides(color = guide_legend(title = NULL))
### This is the DA script:
# Count the number of control and PE cells in the interval
control_count <- sum(subset_data$Condition == conditions[1])
PE_count <- sum(subset_data$Condition == conditions[2])
# Update the density dataframe
density_df$Density_Control[i] <- control_count
density_df$Density_PE[i] <- PE_count
# Count the number of replicates (libraries, samples) per condition
control_samples <- length(unique(subset_data[subset_data$Condition == conditions[1],"sample_id"]))
PE_samples <- length(unique(subset_data[subset_data$Condition == conditions[2],"sample_id"]))
# Update total cells in density dataframe
density_df$n_samples_Control[i] <- control_samples
density_df$n_samples_PE[i] <- PE_samples
# Check if both conditions have enough replications for the test
if (control_samples >= 2 & PE_samples >= 2) {
# Contingency table for Wilcoxon rank-sum test
contingency_table <- table(subset_data$Condition, subset_data$sample_id)
# remove placentas or sn
#contingency_table <- contingency_table[, !grepl("-AB|-sn", colnames(contingency_table))]
# Extract counts for Control and PE conditions
control_v <- as.vector(contingency_table[conditions[1], ])
pe_v <- as.vector(contingency_table[conditions[2], ])
# Remove zeros from the vectors
control_v <- control_v[control_v != 0]
pe_v <- pe_v[pe_v != 0]
# Perform Wilcoxon rank-sum test
test_result <- wilcox.test(pe_v, control_v, exact = F)
# Update Wilcox_Pvalue column in density dataframe
density_df$Wilcox_Pvalue[i] <- test_result$p.value
# Calculate means for PE and Control
Mean_pe <- mean(pe_v)
Mean_control <- mean(control_v)
# Update Sum_Ranks columns in density dataframe
density_df$Mean_Control[i] <- Mean_control
density_df$Mean_PE[i] <- Mean_pe
} else {
# If there are not enough observations for the test, set p-value to 1
density_df$Wilcox_Pvalue[i] <- 1
# Set Sum_Ranks to NA if there are not enough observations
density_df$Mean_PE[i] <- NA
density_df$Mean_Control[i] <- NA
}
}
# Apply Benjamini-Hochberg correction for multiple testing
density_df$Adjusted_Pvalue <- p.adjust(density_df$Wilcox_Pvalue, method = "BH")
# Determine which condition is more abundant in each interval
density_df$More_Abundant <- ifelse(density_df$Density_Control > density_df$Density_PE, conditions[1], conditions[2])
density_df$Difference <- density_df$Density_Control - density_df$Density_PE
# Display significant DA intervals
sig_DA <- density_df[density_df$Adjusted_Pvalue <= 0.2, ]
print(sig_DA)
print(density_df)
# Combine the two variable names with an underscore
if(nrow(sig_DA) >= 1){
new_variable_name <- paste(lineages,gest_age,sep = "_")
sig_DA$test<-new_variable_name
#DA_intervals_df<-rbind(DA_intervals_df,sig_DA)
# Rename the variable
#assign(new_variable_name, sig_DA)
write.csv(sig_DA,paste("~/Documents/ucl_postdoc/cell_trajectories/",lineages,"_",gest_age,".csv", sep = ""))}
# Plot
# Perform correlation tests
cor_expr <- cor.test(density_df$Mean_PE / density_df$Mean_Control, density_df$expr_score)
# Extract correlation and p-value
cor_expr_r <- round(cor_expr$estimate, 2) # Correlation coefficient for IL10
cor_expr_p <- signif(cor_expr$p.value, 2)
# Calculate scaling factor between variables
y1_range <- range(density_df$Mean_PE - density_df$Mean_Control, na.rm = TRUE)
y2_range <- range(density_df$expr_score, na.rm = TRUE)
scaling_factor <- diff(y2_range)/diff(y1_range)
ggplot(data = density_df, aes(x = Pseudotime)) +
geom_vline(xintercept = sig_DA$Pseudotime,
linetype = "solid", color = "lightblue", size = 5) +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +
geom_line(aes(y = Mean_PE - Mean_Control,
color = 'Cells PE - Control'), size = 1) +
geom_line(aes(y = (expr_factor * expr_score)/scaling_factor,
color = 'Expr'), size = 1, linetype = "dotted") +
scale_y_continuous(
name = 'Cells PE - Control',
sec.axis = sec_axis(~ . * scaling_factor, name = "Expression"),
expand = expansion(mult = c(0.1, 0.1)) # Add 10% padding
) +
scale_color_manual(values = c(
'Cells PE - Control' = 'blue',
'Expr' = 'purple'
)) +
labs(
x = "Pseudotime",
title = paste("r =", cor_expr_r, ", p =", cor_expr_p)
) +
theme_minimal() +
theme(legend.position = "top",
plot.title = element_text(size = 10)) +
guides(color = guide_legend(title = NULL))
}Trigger differential abundance analysis
genestoplot <- c("HTRA4","HTRA1","FSTL3","TMEM45A","EGLN3")
seu_placenta1$pseudotime<-seu_placenta1$Lineage1
seu_placenta1$lineages<-"Lineage1"
p1a<-DAbundance_sum_expression_in_pseudotime(seu_placenta1, lineages = "Lineage1", gest_age = NULL,conditions = c("Control", "PE"), pseudotime_intervals = 3,tropho_stb,genestoplot,expr_factor = 1)## No cells found for interval 3 . Skipping...## Pseudotime n_samples_Control n_samples_PE Density_Control Density_PE Wilcox_Pvalue Mean_Control Mean_PE expr_score done Adjusted_Pvalue More_Abundant
## NA NA NA NA NA NA NA NA NA NA <NA> NA <NA>
## 6 15 11 13 514 2149 0.0204785872 46.727273 165.30769 2.199327 yes 0.088740544 PE
## 7 18 11 13 78 924 0.0001997948 7.090909 71.07692 4.186597 yes 0.002597333 PE
## 8 21 10 13 55 467 0.0012288529 5.500000 35.92308 4.583954 yes 0.007987544 PE
## Difference
## NA NA
## 6 -1635
## 7 -846
## 8 -412
## Pseudotime n_samples_Control n_samples_PE Density_Control Density_PE Wilcox_Pvalue Mean_Control Mean_PE expr_score done Adjusted_Pvalue More_Abundant
## 1 0 12 13 228 298 0.6828513737 19.000000 22.92308 0.0867085 yes 1.000000000 PE
## 2 3 12 13 185 618 0.3129806223 15.416667 47.53846 0.2688532 yes 0.993004304 PE
## 3 NA NA NA NA NA NA NA NA NA <NA> NA <NA>
## 4 9 0 2 0 2 1.0000000000 NA NA 2.2951013 yes 1.000000000 PE
## 5 12 12 13 900 2154 0.5137787520 75.000000 165.69231 0.9734493 yes 1.000000000 PE
## 6 15 11 13 514 2149 0.0204785872 46.727273 165.30769 2.1993268 yes 0.088740544 PE
## 7 18 11 13 78 924 0.0001997948 7.090909 71.07692 4.1865974 yes 0.002597333 PE
## 8 21 10 13 55 467 0.0012288529 5.500000 35.92308 4.5839539 yes 0.007987544 PE
## 9 24 12 13 1213 1441 1.0000000000 101.083333 110.84615 2.2257067 yes 1.000000000 PE
## 10 27 12 12 1086 940 0.7287511119 90.500000 78.33333 2.0126593 yes 1.000000000 Control
## 11 30 9 11 1046 740 0.3819247323 116.222222 67.27273 1.5477052 yes 0.993004304 Control
## 12 33 12 11 1488 804 0.8053047818 124.000000 73.09091 0.9222128 yes 1.000000000 Control
## 13 36 12 12 2182 1667 0.8624310005 181.833333 138.91667 0.7114181 yes 1.000000000 Control
## 14 39 11 8 183 103 0.7408525339 16.636364 12.87500 0.3641109 yes 1.000000000 Control
## Difference
## 1 -70
## 2 -433
## 3 NA
## 4 -2
## 5 -1254
## 6 -1635
## 7 -846
## 8 -412
## 9 -228
## 10 146
## 11 306
## 12 684
## 13 515
## 14 80p2a<-DAbundance_sum_expression_in_pseudotime(seu_placenta1, lineages = "Lineage1", conditions = c("Control", "PE"),gest_age = "Early" ,pseudotime_intervals = 3,tropho_stb,genestoplot,expr_factor = 1)## No cells found for interval 3 . Skipping...## Pseudotime n_samples_Control n_samples_PE Density_Control Density_PE Wilcox_Pvalue Mean_Control Mean_PE expr_score done Adjusted_Pvalue More_Abundant
## NA NA NA NA NA NA NA NA NA NA <NA> NA <NA>
## 7 18 4 7 11 513 0.01055554 2.75 73.28571 4.678713 yes 0.1372221 PE
## Difference
## NA NA
## 7 -502
## Pseudotime n_samples_Control n_samples_PE Density_Control Density_PE Wilcox_Pvalue Mean_Control Mean_PE expr_score done Adjusted_Pvalue More_Abundant
## 1 0 5 7 150 74 0.07353002 30.00000 10.571429 0.09639887 yes 0.3186301 Control
## 2 3 5 7 98 173 0.87099119 19.60000 24.714286 0.53818786 yes 0.9435738 PE
## 3 NA NA NA NA NA NA NA NA NA <NA> NA <NA>
## 4 9 0 2 0 2 1.00000000 NA NA 2.29510135 yes 1.0000000 PE
## 5 12 5 7 204 509 0.87099119 40.80000 72.714286 1.08535585 yes 0.9435738 PE
## 6 15 4 7 88 888 0.15637572 22.00000 126.857143 2.76522577 yes 0.5082211 PE
## 7 18 4 7 11 513 0.01055554 2.75000 73.285714 4.67871327 yes 0.1372221 PE
## 8 21 4 7 13 137 0.04621820 3.25000 19.571429 4.80753775 yes 0.3004183 PE
## 9 24 5 7 501 691 0.74490219 100.20000 98.714286 2.17395921 yes 0.9435738 PE
## 10 27 5 6 634 235 0.78225207 126.80000 39.166667 1.54996369 yes 0.9435738 Control
## 11 30 3 5 256 111 0.29090092 85.33333 22.200000 0.90644554 yes 0.6366227 Control
## 12 33 5 5 344 118 0.52580895 68.80000 23.600000 0.64428721 yes 0.8544395 Control
## 13 36 5 6 483 320 0.41024790 96.60000 53.333333 0.48103078 yes 0.7618890 Control
## 14 39 5 3 80 20 0.29382585 16.00000 6.666667 0.31782948 yes 0.6366227 Control
## Difference
## 1 76
## 2 -75
## 3 NA
## 4 -2
## 5 -305
## 6 -800
## 7 -502
## 8 -124
## 9 -190
## 10 399
## 11 145
## 12 226
## 13 163
## 14 60p3a<-DAbundance_sum_expression_in_pseudotime(seu_placenta1, lineages = "Lineage1", conditions = c("Control", "PE"),gest_age = "Late", pseudotime_intervals = 3,tropho_stb,genestoplot,expr_factor = 1)## No cells found for interval 3 . Skipping...## No cells found for interval 4 . Skipping...## Pseudotime n_samples_Control n_samples_PE Density_Control Density_PE Wilcox_Pvalue Mean_Control Mean_PE expr_score done Adjusted_Pvalue More_Abundant
## NA NA NA NA NA NA NA NA NA NA <NA> NA <NA>
## NA.1 NA NA NA NA NA NA NA NA NA <NA> NA <NA>
## 6 15.04979 7 6 421 1257 0.038318763 60.142857 209.50000 1.895712 yes 0.15327505 PE
## 7 18.04979 7 6 67 400 0.011242652 9.571429 66.66667 3.680565 yes 0.06745591 PE
## 8 21.04979 6 6 38 325 0.004998125 6.333333 54.16667 4.526640 yes 0.05997750 PE
## Difference
## NA NA
## NA.1 NA
## 6 -836
## 7 -333
## 8 -287
## Pseudotime n_samples_Control n_samples_PE Density_Control Density_PE Wilcox_Pvalue Mean_Control Mean_PE expr_score done Adjusted_Pvalue More_Abundant
## 1 0.0497923 7 6 79 230 0.174143130 11.285714 38.33333 0.08234964 yes 0.41794351 PE
## 2 3.0497923 7 6 86 439 0.429488462 12.285714 73.16667 0.13068479 yes 0.73626593 PE
## 3 NA NA NA NA NA NA NA NA NA <NA> NA <NA>
## 4 NA NA NA NA NA NA NA NA NA <NA> NA <NA>
## 5 12.0497923 7 6 705 1665 0.224638639 100.714286 277.50000 0.94051534 yes 0.44927728 PE
## 6 15.0497923 7 6 421 1257 0.038318763 60.142857 209.50000 1.89571236 yes 0.15327505 PE
## 7 18.0497923 7 6 67 400 0.011242652 9.571429 66.66667 3.68056458 yes 0.06745591 PE
## 8 21.0497923 6 6 38 325 0.004998125 6.333333 54.16667 4.52663981 yes 0.05997750 PE
## 9 24.0497923 7 6 718 758 0.617075077 102.571429 126.33333 2.26724496 yes 0.86474027 PE
## 10 27.0497923 7 6 457 707 0.133614403 65.285714 117.83333 2.36135889 yes 0.40084321 PE
## 11 30.0497923 6 6 794 626 0.936186293 132.333333 104.33333 1.69935966 yes 0.93618629 Control
## 12 33.0497923 7 6 1166 706 0.720616893 166.571429 117.66667 0.98173188 yes 0.86474027 Control
## 13 36.0497923 7 6 1678 1331 0.830324258 239.714286 221.83333 0.76887975 yes 0.90580828 Control
## 14 39.0497923 6 5 87 72 0.713753631 14.500000 14.40000 0.36973994 yes 0.86474027 Control
## Difference
## 1 -151
## 2 -353
## 3 NA
## 4 NA
## 5 -960
## 6 -836
## 7 -333
## 8 -287
## 9 -40
## 10 -250
## 11 168
## 12 460
## 13 347
## 14 15p1a+p2a+p3a## Warning: Removed 1 row containing missing values or values outside the scale range (`geom_vline()`).## Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).
## Removed 1 row containing missing values or values outside the scale range (`geom_line()`).## Warning: Removed 1 row containing missing values or values outside the scale range (`geom_vline()`).## Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).
## Removed 1 row containing missing values or values outside the scale range (`geom_line()`).## Warning: Removed 2 rows containing missing values or values outside the scale range (`geom_vline()`).## Warning: Removed 2 rows containing missing values or values outside the scale range (`geom_line()`).
## Removed 2 rows containing missing values or values outside the scale range (`geom_line()`).
p2a+p3a## Warning: Removed 1 row containing missing values or values outside the scale range (`geom_vline()`).## Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).
## Removed 1 row containing missing values or values outside the scale range (`geom_line()`).## Warning: Removed 2 rows containing missing values or values outside the scale range (`geom_vline()`).## Warning: Removed 2 rows containing missing values or values outside the scale range (`geom_line()`).
## Removed 2 rows containing missing values or values outside the scale range (`geom_line()`).
Function to calculate differential abundance in proportion to total
trophoblasts per donor
Proportion_DAbundance_sum_expression_in_pseudotime <- function(df_sample, lineages = "Lineage1", gest_age = NULL, conditions = c("Control", "PE"), pseudotime_intervals = 1, tropho,genestoplot,expr_factor = 100) {
# Subset data frame by lineages and gestational age if specified
### make the script:
#df_sample<-seu_placenta1
#lineages = "Lineage1"
#gest_age = NULL
#conditions = c("Control", "PE")
#pseudotime_intervals = 3
#tropho<-tropho_stb
#expr_factor = 25
total_trophos_sample<-table(df_sample$sample_id)
total_trophos_sample_df <- as.data.frame(t(as.matrix(total_trophos_sample)))
if (!is.null(gest_age)) {
dfl1 <- df_sample[df_sample$lineages == lineages & df_sample$GA_Category == gest_age, ]
} else {
dfl1 <- df_sample[df_sample$lineages == lineages, ]
}
# Define the range of pseudotime values
min_pt <- min(dfl1$pseudotime)
max_pt <- max(dfl1$pseudotime)
# Create a sequence of pseudotime values at intervals of specified units
pseudotime_seq <- seq(min_pt, max_pt, by = pseudotime_intervals)
# Initialize density_df outside the loop if not already done
# Initialize empty dataframe to store density values for control and PE samples
density_df <- data.frame(Pseudotime = pseudotime_seq,
n_samples_Control = rep(0, length(pseudotime_seq)),
n_samples_PE = rep(0, length(pseudotime_seq)),
Density_Control = rep(0, length(pseudotime_seq)),
Density_PE = rep(0, length(pseudotime_seq)),
Wilcox_Pvalue = rep(0, length(pseudotime_seq)),
Mean_Control = rep(0, length(pseudotime_seq)),
Mean_PE = rep(0, length(pseudotime_seq)),
expr_score = NA,
done = "no"
)
# Loop through pseudotime sequence
for (i in seq_along(pseudotime_seq)) {
# Subset data for each pseudotime interval
subset_data <- dfl1[dfl1$pseudotime >= pseudotime_seq[i] & dfl1$pseudotime < pseudotime_seq[i] + pseudotime_intervals, ]
# Check if there are any cells in this interval
if (nrow(subset_data) == 0 || length(intersect(rownames(subset_data), colnames(tropho))) == 0) {
message(paste("No cells found for interval", i, ". Skipping..."))
density_df$done[i] <- "yes"
density_df[i,]<-NA
next # Skip to the next iteration
}
# Subset the Seurat object
interval_seurat <- subset(tropho, cells = intersect(rownames(subset_data), colnames(tropho)))
# Check if the subset Seurat object has any cells
if (ncol(interval_seurat) <= 1) {
message(paste("No cells remained after subsetting for interval", i, ". Skipping..."))
next # Skip to the next iteration
}
#genes_of_interest <- list(c("LEP", "LTF", "DERL3"))
#interval_seurat <- AddModuleScore(interval_seurat,
# features = genes_of_interest,
# name = "GeneScore")
#density_df$expr_score[i] <- mean(interval_seurat$GeneScore1)
# Get the expression data for genes of interest
#genes_of_interest <- c("HTRA1","FSTL3","LTF","LEP","FLT1")
genes_of_interest <- genestoplot
gene_expr <- GetAssayData(interval_seurat, slot = "data")[genes_of_interest, ]
# Calculate the sum of expression of all the genes for each cell
if(length(genes_of_interest) == 1 ) {
sum_expr <- t(as.data.frame(gene_expr))} else {
sum_expr <- colSums(gene_expr)}
# Calculate the mean of the sum across all cells in this interval
density_df$expr_score[i] <- mean(sum_expr)
density_df$done[i] <- "yes"
#}
# ggplot(data = density_df, aes(x = pseudotime)) +
# geom_line(aes(y = 10*expr_score, color = 'coral'), size = 1) +
# labs(x = "Pseudotime", y = "Mean Hypoxia score", title = paste(lineages,"and",gest_age)) +
# theme_minimal() +
# theme(legend.position = "top") +
# guides(color = guide_legend(title = NULL))
### This is the DA script:
# Count the number of control and PE cells in the interval
control_count <- sum(subset_data$Condition == conditions[1])
PE_count <- sum(subset_data$Condition == conditions[2])
# Update the density dataframe
density_df$Density_Control[i] <- control_count
density_df$Density_PE[i] <- PE_count
# Count the number of replicates (libraries, samples) per condition
control_samples <- length(unique(subset_data[subset_data$Condition == conditions[1],"sample_id"]))
PE_samples <- length(unique(subset_data[subset_data$Condition == conditions[2],"sample_id"]))
# Update total cells in density dataframe
density_df$n_samples_Control[i] <- control_samples
density_df$n_samples_PE[i] <- PE_samples
# Check if both conditions have enough replications for the test
if (control_samples >= 2 & PE_samples >= 2) {
# Contingency table for Wilcoxon rank-sum test
contingency_table <- table(subset_data$Condition, subset_data$sample_id)
# remove placentas or sn
#contingency_table <- contingency_table[, !grepl("-AB|-sn", colnames(contingency_table))]
contingency_df <- as.data.frame.matrix(contingency_table)
# Ensure only the common columns are kept
common_cols <- intersect(colnames(contingency_df), colnames(total_trophos_sample_df))
total_trophos_sample_filtered <- total_trophos_sample_df %>% select(all_of(common_cols))
# Merge by row binding
merged_table <- rbind(contingency_df, total_trophos_sample_filtered)
merged_table["C_prop",]<-merged_table[1,]/merged_table[3,]
merged_table["PE_prop",]<-merged_table[2,]/merged_table[3,]
# Extract counts for Control and PE conditions
control_v <- as.vector(unlist(merged_table["C_prop",]))
pe_v <- as.vector(unlist(merged_table["PE_prop",]))
# Remove zeros from the vectors
control_v <- control_v[control_v != 0]
pe_v <- pe_v[pe_v != 0]
# Perform Wilcoxon rank-sum test
test_result <- wilcox.test(pe_v, control_v, exact = F)
# Update Wilcox_Pvalue column in density dataframe
density_df$Wilcox_Pvalue[i] <- test_result$p.value
# Calculate means for PE and Control
Mean_pe <- mean(pe_v)
Mean_control <- mean(control_v)
# Update Sum_Ranks columns in density dataframe
density_df$Mean_Control[i] <- Mean_control
density_df$Mean_PE[i] <- Mean_pe
} else {
# If there are not enough observations for the test, set p-value to 1
density_df$Wilcox_Pvalue[i] <- 1
# Set Sum_Ranks to NA if there are not enough observations
density_df$Mean_PE[i] <- NA
density_df$Mean_Control[i] <- NA
}
}
# Apply Benjamini-Hochberg correction for multiple testing
density_df$Adjusted_Pvalue <- p.adjust(density_df$Wilcox_Pvalue, method = "BH")
# Determine which condition is more abundant in each interval
density_df$More_Abundant <- ifelse(density_df$Density_Control > density_df$Density_PE, conditions[1], conditions[2])
density_df$Difference <- density_df$Density_Control - density_df$Density_PE
# Display significant DA intervals
sig_DA <- density_df[density_df$Adjusted_Pvalue <= 0.2, ]
print(sig_DA)
print(density_df)
# Combine the two variable names with an underscore
if(nrow(sig_DA) >= 1){
new_variable_name <- paste(lineages,gest_age,sep = "_")
sig_DA$test<-new_variable_name
#DA_intervals_df<-rbind(DA_intervals_df,sig_DA)
# Rename the variable
#assign(new_variable_name, sig_DA)
write.csv(sig_DA,paste("~/Documents/ucl_postdoc/cell_trajectories/",lineages,"_",gest_age,".csv", sep = ""))}
# Plot
# Perform correlation tests
cor_expr <- cor.test(density_df$Mean_PE / density_df$Mean_Control, density_df$expr_score)
# Extract correlation and p-value
cor_expr_r <- round(cor_expr$estimate, 2) # Correlation coefficient for IL10
cor_expr_p <- signif(cor_expr$p.value, 2)
#expr_factor = 1
# Add correlation results to the plot title
# Calculate scaling factor between variables
y1_range <- range(density_df$Mean_PE - density_df$Mean_Control, na.rm = TRUE)
y2_range <- range(density_df$expr_score, na.rm = TRUE)
scaling_factor <- diff(y2_range)/diff(y1_range)
ggplot(data = density_df, aes(x = Pseudotime)) +
geom_vline(xintercept = sig_DA$Pseudotime,
linetype = "solid", color = "lightblue", size = 5) +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +
geom_line(aes(y = Mean_PE - Mean_Control,
color = 'Proportion Cells PE - Control'), size = 1) +
geom_line(aes(y = (expr_score)/scaling_factor,
color = 'Expr'), size = 1, linetype = "dotted") +
scale_y_continuous(
name = 'Proportion Cells PE - Control',
sec.axis = sec_axis(~ . * scaling_factor, name = "Expression"),
expand = expansion(mult = c(0.1, 0.1)) # Add 10% padding
) +
scale_color_manual(values = c(
'Proportion Cells PE - Control' = 'blue',
'Expr' = 'purple'
)) +
labs(
x = "Pseudotime",
title = paste("r =", cor_expr_r, ", p =", cor_expr_p)
) +
theme_minimal() +
theme(legend.position = "top",
plot.title = element_text(size = 10)) +
guides(color = guide_legend(title = NULL))
}genestoplot <- c("HTRA4","HTRA1","FSTL3","TMEM45A","EGLN3")
seu_placenta1$pseudotime<-seu_placenta1$Lineage1
seu_placenta1$lineages<-"Lineage1"
p1a<-Proportion_DAbundance_sum_expression_in_pseudotime(seu_placenta1, lineages = "Lineage1", gest_age = NULL,conditions = c("Control", "PE"), pseudotime_intervals = 3,tropho_stb,genestoplot,expr_factor = 1)## No cells found for interval 3 . Skipping...## Pseudotime n_samples_Control n_samples_PE Density_Control Density_PE Wilcox_Pvalue Mean_Control Mean_PE expr_score done Adjusted_Pvalue More_Abundant
## NA NA NA NA NA NA NA NA NA NA <NA> NA <NA>
## 6 15 11 13 514 2149 0.0021360271 0.06711436 0.178501573 2.1993268 yes 0.009256118 PE
## 7 18 11 13 78 924 0.0004091209 0.01457911 0.093995155 4.1865974 yes 0.002979325 PE
## 8 21 10 13 55 467 0.0004583576 0.01080712 0.042839312 4.5839539 yes 0.002979325 PE
## 12 33 12 11 1488 804 0.0604981907 0.12416864 0.057211526 0.9222128 yes 0.157295296 Control
## 14 39 11 8 183 103 0.0430709909 0.02443279 0.008779221 0.3641109 yes 0.139980720 Control
## Difference
## NA NA
## 6 -1635
## 7 -846
## 8 -412
## 12 684
## 14 80
## Pseudotime n_samples_Control n_samples_PE Density_Control Density_PE Wilcox_Pvalue Mean_Control Mean_PE expr_score done Adjusted_Pvalue More_Abundant
## 1 0 12 13 228 298 0.4627639706 0.08000466 0.025715500 0.0867085 yes 0.751991452 PE
## 2 3 12 13 185 618 0.6833131587 0.04793879 0.125269165 0.2688532 yes 0.807551915 PE
## 3 NA NA NA NA NA NA NA NA NA <NA> NA <NA>
## 4 9 0 2 0 2 1.0000000000 NA NA 2.2951013 yes 1.000000000 PE
## 5 12 12 13 900 2154 0.6438382008 0.10819060 0.129218095 0.9734493 yes 0.807551915 PE
## 6 15 11 13 514 2149 0.0021360271 0.06711436 0.178501573 2.1993268 yes 0.009256118 PE
## 7 18 11 13 78 924 0.0004091209 0.01457911 0.093995155 4.1865974 yes 0.002979325 PE
## 8 21 10 13 55 467 0.0004583576 0.01080712 0.042839312 4.5839539 yes 0.002979325 PE
## 9 24 12 13 1213 1441 0.6438382008 0.17306431 0.121967513 2.2257067 yes 0.807551915 PE
## 10 27 12 12 1086 940 0.8852339145 0.11221777 0.078759310 2.0126593 yes 0.959003407 Control
## 11 30 9 11 1046 740 0.1965119202 0.10510086 0.062430217 1.5477052 yes 0.364950709 Control
## 12 33 12 11 1488 804 0.0604981907 0.12416864 0.057211526 0.9222128 yes 0.157295296 Control
## 13 36 12 12 2182 1667 0.1572127533 0.16930124 0.109778096 0.7114181 yes 0.340627632 Control
## 14 39 11 8 183 103 0.0430709909 0.02443279 0.008779221 0.3641109 yes 0.139980720 Control
## Difference
## 1 -70
## 2 -433
## 3 NA
## 4 -2
## 5 -1254
## 6 -1635
## 7 -846
## 8 -412
## 9 -228
## 10 146
## 11 306
## 12 684
## 13 515
## 14 80p2a<-Proportion_DAbundance_sum_expression_in_pseudotime(seu_placenta1, lineages = "Lineage1", conditions = c("Control", "PE"),gest_age = "Early" ,pseudotime_intervals = 3,tropho_stb,genestoplot,expr_factor = 1)## No cells found for interval 3 . Skipping...## Pseudotime n_samples_Control n_samples_PE Density_Control Density_PE Wilcox_Pvalue Mean_Control Mean_PE expr_score done Adjusted_Pvalue More_Abundant
## NA NA NA NA NA NA NA NA NA NA <NA> NA <NA>
## 6 15 4 7 88 888 0.01073342 0.029063160 0.213684048 2.7652258 yes 0.06976720 PE
## 7 18 4 7 11 513 0.01073342 0.005105015 0.134150261 4.6787133 yes 0.06976720 PE
## 8 21 4 7 13 137 0.01816302 0.007364155 0.044325094 4.8075377 yes 0.07870641 PE
## 14 39 5 3 80 20 0.03688843 0.035400813 0.007374697 0.3178295 yes 0.11988738 Control
## Difference
## NA NA
## 6 -800
## 7 -502
## 8 -124
## 14 60
## Pseudotime n_samples_Control n_samples_PE Density_Control Density_PE Wilcox_Pvalue Mean_Control Mean_PE expr_score done Adjusted_Pvalue More_Abundant
## 1 0 5 7 150 74 0.25562311 0.159818261 0.025416608 0.09639887 yes 0.41538755 Control
## 2 3 5 7 98 173 0.74533307 0.082733130 0.188390416 0.53818786 yes 0.96893299 PE
## 3 NA NA NA NA NA NA NA NA NA <NA> NA <NA>
## 4 9 0 2 0 2 1.00000000 NA NA 2.29510135 yes 1.00000000 PE
## 5 12 5 7 204 509 0.62611748 0.091496164 0.097974205 1.08535585 yes 0.90439192 PE
## 6 15 4 7 88 888 0.01073342 0.029063160 0.213684048 2.76522577 yes 0.06976720 PE
## 7 18 4 7 11 513 0.01073342 0.005105015 0.134150261 4.67871327 yes 0.06976720 PE
## 8 21 4 7 13 137 0.01816302 0.007364155 0.044325094 4.80753775 yes 0.07870641 PE
## 9 24 5 7 501 691 0.87099119 0.184983211 0.146129228 2.17395921 yes 1.00000000 PE
## 10 27 5 6 634 235 1.00000000 0.120461053 0.045396145 1.54996369 yes 1.00000000 Control
## 11 30 3 5 256 111 0.23303798 0.079148873 0.021416543 0.90644554 yes 0.41538755 Control
## 12 33 5 5 344 118 0.21007504 0.096996469 0.031627849 0.64428721 yes 0.41538755 Control
## 13 36 5 6 483 320 0.08283743 0.147395710 0.077685325 0.48103078 yes 0.21537731 Control
## 14 39 5 3 80 20 0.03688843 0.035400813 0.007374697 0.31782948 yes 0.11988738 Control
## Difference
## 1 76
## 2 -75
## 3 NA
## 4 -2
## 5 -305
## 6 -800
## 7 -502
## 8 -124
## 9 -190
## 10 399
## 11 145
## 12 226
## 13 163
## 14 60p3a<-Proportion_DAbundance_sum_expression_in_pseudotime(seu_placenta1, lineages = "Lineage1", conditions = c("Control", "PE"),gest_age = "Late", pseudotime_intervals = 3,tropho_stb,genestoplot,expr_factor = 1)## No cells found for interval 3 . Skipping...## No cells found for interval 4 . Skipping...## Pseudotime n_samples_Control n_samples_PE Density_Control Density_PE Wilcox_Pvalue Mean_Control Mean_PE expr_score done Adjusted_Pvalue More_Abundant
## NA NA NA NA NA NA NA NA NA NA <NA> NA <NA>
## NA.1 NA NA NA NA NA NA NA NA NA <NA> NA <NA>
## 8 21.04979 6 6 38 325 0.01306523 0.0117548 0.04059447 4.52664 yes 0.1567827 PE
## Difference
## NA NA
## NA.1 NA
## 8 -287
## Pseudotime n_samples_Control n_samples_PE Density_Control Density_PE Wilcox_Pvalue Mean_Control Mean_PE expr_score done Adjusted_Pvalue More_Abundant
## 1 0.0497923 7 6 79 230 0.83032426 0.02316665 0.026688460 0.08234964 yes 0.9058083 PE
## 2 3.0497923 7 6 86 439 0.83032426 0.02291398 0.051003454 0.13068479 yes 0.9058083 PE
## 3 NA NA NA NA NA NA NA NA NA <NA> NA <NA>
## 4 NA NA NA NA NA NA NA NA NA <NA> NA <NA>
## 5 12.0497923 7 6 705 1665 0.83032426 0.12152602 0.167894545 0.94051534 yes 0.9058083 PE
## 6 15.0497923 7 6 421 1257 0.10041249 0.08886009 0.137223242 1.89571236 yes 0.4016500 PE
## 7 18.0497923 7 6 67 400 0.03831876 0.01973499 0.045665823 3.68056458 yes 0.2299126 PE
## 8 21.0497923 6 6 38 325 0.01306523 0.01175480 0.040594473 4.52663981 yes 0.1567827 PE
## 9 24.0497923 7 6 718 758 0.72098486 0.16593457 0.094771129 2.26724496 yes 0.9058083 PE
## 10 27.0497923 7 6 457 707 0.94305667 0.10635471 0.113021014 2.36135889 yes 0.9430567 PE
## 11 30.0497923 6 6 794 626 0.81018124 0.11952059 0.095424715 1.69935966 yes 0.9058083 Control
## 12 33.0497923 7 6 1166 706 0.17473582 0.14573872 0.080810465 0.98173188 yes 0.5242075 Control
## 13 36.0497923 7 6 1678 1331 0.83032426 0.18196133 0.140155744 0.76887975 yes 0.9058083 Control
## 14 39.0497923 6 5 87 72 0.52281665 0.01316838 0.008096325 0.36973994 yes 0.9058083 Control
## Difference
## 1 -151
## 2 -353
## 3 NA
## 4 NA
## 5 -960
## 6 -836
## 7 -333
## 8 -287
## 9 -40
## 10 -250
## 11 168
## 12 460
## 13 347
## 14 15p2a+p3a## Warning: Removed 1 row containing missing values or values outside the scale range (`geom_vline()`).## Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).
## Removed 1 row containing missing values or values outside the scale range (`geom_line()`).## Warning: Removed 2 rows containing missing values or values outside the scale range (`geom_vline()`).## Warning: Removed 2 rows containing missing values or values outside the scale range (`geom_line()`).
## Removed 2 rows containing missing values or values outside the scale range (`geom_line()`).
For main figure:
seu_placenta1
seu_placenta2<-subset(seu_placenta1, umap_1 < 5)
# Create the initial ggplot
p <- ggplot(seu_placenta2, aes(x = umap_1, y = umap_2, col = scaled_scores)) +
geom_point() +
scale_color_viridis_c(option = "C")
# Add the geom_path layer to the ggplot
curves <- slingCurves(seu_placenta, as.df = TRUE)
curves2<-subset(curves, Order > 38)
p1 <- p +
geom_path(data = curves2 %>% arrange(Order),
aes(group = Lineage), col = "green", size = 1.5,
arrow = arrow(type = "closed", length = unit(0.15, "inches"))) +
theme_classic() +
labs(col = "Lineage") +
#scale_color_brewer(palette = "Set1")+
geom_point(data = curves2 %>%
group_by(Lineage) %>%
slice(round(seq(1, n(), length.out = 7))),
aes( size = 4),col = "green") +
theme_classic() +
labs(col = "Lineage")
#scale_color_brewer(palette = "Set1")
p1
#ppp4<-ggplot(seu_placenta1, aes(x = umap_1, y = umap_2, col = scaled_scores)) +
# geom_point() +
# scale_color_viridis_c(option = "C")
#ppp4genestoplot <- c("HTRA4","HTRA1","FSTL3","TMEM45A","EGLN3")
seu_placenta2$pseudotime<-seu_placenta2$Lineage1
seu_placenta2$lineages<-"Lineage1"
p1a<-DAbundance_sum_expression_in_pseudotime(seu_placenta2, lineages = "Lineage1", gest_age = NULL,conditions = c("Control", "PE"), pseudotime_intervals = 3,tropho_stb,genestoplot,expr_factor = 1)## Pseudotime n_samples_Control n_samples_PE Density_Control Density_PE Wilcox_Pvalue Mean_Control Mean_PE expr_score done Adjusted_Pvalue More_Abundant
## 2 14.7469 11 13 599 2287 0.0276943716 54.454545 175.92308 2.033027 yes 0.092314572 PE
## 3 17.7469 11 13 91 1007 0.0003801491 8.272727 77.46154 4.042130 yes 0.003801491 PE
## 4 20.7469 10 13 62 515 0.0010798255 6.200000 39.61538 4.533899 yes 0.005399128 PE
## Difference
## 2 -1688
## 3 -916
## 4 -453
## Pseudotime n_samples_Control n_samples_PE Density_Control Density_PE Wilcox_Pvalue Mean_Control Mean_PE expr_score done Adjusted_Pvalue More_Abundant
## 1 11.7469 12 13 795 1886 0.4462712558 66.250000 145.07692 0.9307012 yes 0.763925505 PE
## 2 14.7469 11 13 599 2287 0.0276943716 54.454545 175.92308 2.0330270 yes 0.092314572 PE
## 3 17.7469 11 13 91 1007 0.0003801491 8.272727 77.46154 4.0421301 yes 0.003801491 PE
## 4 20.7469 10 13 62 515 0.0010798255 6.200000 39.61538 4.5338991 yes 0.005399128 PE
## 5 23.7469 12 13 1095 1379 0.8065268415 91.250000 106.07692 2.2523024 yes 0.896140935 PE
## 6 26.7469 12 12 1142 912 0.7505170643 95.166667 76.00000 1.9937810 yes 0.896140935 Control
## 7 29.7469 9 11 1031 796 0.4472408649 114.555556 72.36364 1.6088530 yes 0.763925505 Control
## 8 32.7469 12 11 1352 697 0.5587598229 112.666667 63.36364 0.9525701 yes 0.798228318 Control
## 9 35.7469 12 12 2268 1709 0.9080327383 189.000000 142.41667 0.7414424 yes 0.908032738 Control
## 10 38.7469 11 10 310 203 0.4583553030 28.181818 20.30000 0.3955224 yes 0.763925505 Control
## Difference
## 1 -1091
## 2 -1688
## 3 -916
## 4 -453
## 5 -284
## 6 230
## 7 235
## 8 655
## 9 559
## 10 107p2a<-DAbundance_sum_expression_in_pseudotime(seu_placenta2, lineages = "Lineage1", conditions = c("Control", "PE"),gest_age = "Early" ,pseudotime_intervals = 3,tropho_stb,genestoplot,expr_factor = 1)## Pseudotime n_samples_Control n_samples_PE Density_Control Density_PE Wilcox_Pvalue Mean_Control Mean_PE expr_score done Adjusted_Pvalue More_Abundant
## 3 17.7469 4 7 12 555 0.01358654 3.0 79.28571 4.603160 yes 0.1358654 PE
## 4 20.7469 4 7 14 155 0.03547995 3.5 22.14286 4.718515 yes 0.1773997 PE
## Difference
## 3 -543
## 4 -141
## Pseudotime n_samples_Control n_samples_PE Density_Control Density_PE Wilcox_Pvalue Mean_Control Mean_PE expr_score done Adjusted_Pvalue More_Abundant
## 1 11.7469 5 7 187 442 0.87099119 37.40000 63.14286 0.9911097 yes 0.9677680 PE
## 2 14.7469 4 7 103 897 0.15637572 25.75000 128.14286 2.6192571 yes 0.3909393 PE
## 3 17.7469 4 7 12 555 0.01358654 3.00000 79.28571 4.6031599 yes 0.1358654 PE
## 4 20.7469 4 7 14 155 0.03547995 3.50000 22.14286 4.7185146 yes 0.1773997 PE
## 5 23.7469 5 7 434 672 1.00000000 86.80000 96.00000 2.2336472 yes 1.0000000 PE
## 6 26.7469 5 6 678 240 0.78225207 135.60000 40.00000 1.5439314 yes 0.9677680 Control
## 7 29.7469 3 5 259 119 0.54859547 86.33333 23.80000 0.9766815 yes 0.7837078 Control
## 8 32.7469 5 5 326 95 0.40339530 65.20000 19.00000 0.5942921 yes 0.6837465 Control
## 9 35.7469 5 6 478 337 0.41024790 95.60000 56.16667 0.5388545 yes 0.6837465 Control
## 10 38.7469 5 4 123 32 0.11134689 24.60000 8.00000 0.2604051 yes 0.3711563 Control
## Difference
## 1 -255
## 2 -794
## 3 -543
## 4 -141
## 5 -238
## 6 438
## 7 140
## 8 231
## 9 141
## 10 91p3a<-DAbundance_sum_expression_in_pseudotime(seu_placenta2, lineages = "Lineage1", conditions = c("Control", "PE"),gest_age = "Late", pseudotime_intervals = 3,tropho_stb,genestoplot,expr_factor = 1)## Pseudotime n_samples_Control n_samples_PE Density_Control Density_PE Wilcox_Pvalue Mean_Control Mean_PE expr_score done Adjusted_Pvalue More_Abundant
## 2 15.88236 7 6 277 1009 0.012299213 39.571429 168.16667 2.36611 yes 0.05534646 PE
## 3 18.88236 7 6 50 383 0.006273708 7.142857 63.83333 4.03582 yes 0.05534646 PE
## Difference
## 2 -732
## 3 -333
## Pseudotime n_samples_Control n_samples_PE Density_Control Density_PE Wilcox_Pvalue Mean_Control Mean_PE expr_score done Adjusted_Pvalue More_Abundant
## 1 12.88236 7 6 885 2076 0.224638639 126.428571 346.00000 1.0350638 yes 0.37866788 PE
## 2 15.88236 7 6 277 1009 0.012299213 39.571429 168.16667 2.3661098 yes 0.05534646 PE
## 3 18.88236 7 6 50 383 0.006273708 7.142857 63.83333 4.0358196 yes 0.05534646 PE
## 4 21.88236 7 6 595 815 0.252445252 85.000000 135.83333 2.6044763 yes 0.37866788 PE
## 5 24.88236 7 6 291 288 0.174735823 41.571429 48.00000 2.2077441 yes 0.37866788 Control
## 6 27.88236 7 6 538 793 0.224638639 76.857143 132.16667 2.2929353 yes 0.37866788 PE
## 7 30.88236 7 6 743 474 0.617075077 106.142857 79.00000 1.4327672 yes 0.69420946 Control
## 8 33.88236 7 6 1469 1043 0.520316801 209.857143 173.83333 0.9702918 yes 0.66897874 Control
## 9 36.88236 7 6 1283 966 0.943056671 183.285714 161.00000 0.6454614 yes 0.94305667 Control
## Difference
## 1 -1191
## 2 -732
## 3 -333
## 4 -220
## 5 3
## 6 -255
## 7 269
## 8 426
## 9 317p1a+p2a+p3a
p2a+p3a
genestoplot <- c("HTRA4","HTRA1","FSTL3","TMEM45A","EGLN3")
seu_placenta1$pseudotime<-seu_placenta1$Lineage1
seu_placenta1$lineages<-"Lineage1"
#p1a<-Proportion_DAbundance_sum_expression_in_pseudotime(seu_placenta2, lineages = "Lineage1", gest_age = NULL,conditions = c("Control", "PE"), pseudotime_intervals = 3,tropho_stb,genestoplot,expr_factor = 0.1)
p2a<-Proportion_DAbundance_sum_expression_in_pseudotime(seu_placenta2, lineages = "Lineage1", conditions = c("Control", "PE"),gest_age = "Early" ,pseudotime_intervals = 3,tropho_stb,genestoplot,expr_factor = 1)## Pseudotime n_samples_Control n_samples_PE Density_Control Density_PE Wilcox_Pvalue Mean_Control Mean_PE expr_score done Adjusted_Pvalue More_Abundant
## 2 14.7469 4 7 103 897 0.01073342 0.040785302 0.277044432 2.6192571 yes 0.05366708 PE
## 3 17.7469 4 7 12 555 0.01073342 0.006799318 0.192204471 4.6031599 yes 0.05366708 PE
## 4 20.7469 4 7 14 155 0.02975807 0.012883839 0.083213953 4.7185146 yes 0.07342771 PE
## 8 32.7469 5 5 326 95 0.03671386 0.102247716 0.024603369 0.5942921 yes 0.07342771 Control
## 9 35.7469 5 6 478 337 0.05523425 0.200614206 0.095327414 0.5388545 yes 0.09205709 Control
## 10 38.7469 5 4 123 32 0.01996445 0.081029512 0.008507792 0.2604051 yes 0.06654818 Control
## Difference
## 2 -794
## 3 -543
## 4 -141
## 8 231
## 9 141
## 10 91
## Pseudotime n_samples_Control n_samples_PE Density_Control Density_PE Wilcox_Pvalue Mean_Control Mean_PE expr_score done Adjusted_Pvalue More_Abundant
## 1 11.7469 5 7 187 442 1.00000000 0.142330078 0.110557671 0.9911097 yes 1.00000000 PE
## 2 14.7469 4 7 103 897 0.01073342 0.040785302 0.277044432 2.6192571 yes 0.05366708 PE
## 3 17.7469 4 7 12 555 0.01073342 0.006799318 0.192204471 4.6031599 yes 0.05366708 PE
## 4 20.7469 4 7 14 155 0.02975807 0.012883839 0.083213953 4.7185146 yes 0.07342771 PE
## 5 23.7469 5 7 434 672 0.87099119 0.239458971 0.167855265 2.2336472 yes 1.00000000 PE
## 6 26.7469 5 6 678 240 0.92726447 0.137127312 0.055494053 1.5439314 yes 1.00000000 Control
## 7 29.7469 3 5 259 119 0.55098499 0.081362396 0.024378527 0.9766815 yes 0.78712141 Control
## 8 32.7469 5 5 326 95 0.03671386 0.102247716 0.024603369 0.5942921 yes 0.07342771 Control
## 9 35.7469 5 6 478 337 0.05523425 0.200614206 0.095327414 0.5388545 yes 0.09205709 Control
## 10 38.7469 5 4 123 32 0.01996445 0.081029512 0.008507792 0.2604051 yes 0.06654818 Control
## Difference
## 1 -255
## 2 -794
## 3 -543
## 4 -141
## 5 -238
## 6 438
## 7 140
## 8 231
## 9 141
## 10 91p3a<-Proportion_DAbundance_sum_expression_in_pseudotime(seu_placenta2, lineages = "Lineage1", conditions = c("Control", "PE"),gest_age = "Late", pseudotime_intervals = 3,tropho_stb,genestoplot,expr_factor = 1)## Pseudotime n_samples_Control n_samples_PE Density_Control Density_PE Wilcox_Pvalue Mean_Control Mean_PE expr_score done Adjusted_Pvalue More_Abundant
## 2 15.88236 7 6 277 1009 0.03831876 0.07238716 0.12482628 2.36611 yes 0.1724344 PE
## 3 18.88236 7 6 50 383 0.02680913 0.01662058 0.04841309 4.03582 yes 0.1724344 PE
## Difference
## 2 -732
## 3 -333
## Pseudotime n_samples_Control n_samples_PE Density_Control Density_PE Wilcox_Pvalue Mean_Control Mean_PE expr_score done Adjusted_Pvalue More_Abundant
## 1 12.88236 7 6 885 2076 0.72098486 0.16179616 0.24251174 1.0350638 yes 0.8111080 PE
## 2 15.88236 7 6 277 1009 0.03831876 0.07238716 0.12482628 2.3661098 yes 0.1724344 PE
## 3 18.88236 7 6 50 383 0.02680913 0.01662058 0.04841309 4.0358196 yes 0.1724344 PE
## 4 21.88236 7 6 595 815 0.83032426 0.13813626 0.10623763 2.6044763 yes 0.8303243 PE
## 5 24.88236 7 6 291 288 0.72098486 0.08711023 0.03922031 2.2077441 yes 0.8111080 Control
## 6 27.88236 7 6 538 793 0.72098486 0.09819989 0.14112687 2.2929353 yes 0.8111080 PE
## 7 30.88236 7 6 743 474 0.52031680 0.10032333 0.07009519 1.4327672 yes 0.8111080 Control
## 8 33.88236 7 6 1469 1043 0.35311123 0.17707654 0.12417391 0.9702918 yes 0.8111080 Control
## 9 36.88236 7 6 1283 966 0.43203489 0.14834985 0.10339498 0.6454614 yes 0.8111080 Control
## Difference
## 1 -1191
## 2 -732
## 3 -333
## 4 -220
## 5 3
## 6 -255
## 7 269
## 8 426
## 9 317#p1a+p2a+p3a
p2a+p3a